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A  three  node  flat  shell  element  with  six  engineering  displacement 
degrees-o£-f reedom  at  each  node  is  developed.  The  basic 
formulation  allows  for  the  arbitrary  location  of  the  reference 
surface  in  which  the  membrane  forces  and  bending  moments  are  fully 
coupled. 

The  well-known,  highly  accurate,  OKT  bending  element  is  combined 
with  a  higher  order  membrane  element  in  order  to  obtain  a 
consistent  formulation.  The  higher  order  membrane  behavior  is 
obtained  by  the  introduction  of  three  additional  normal  rotational 
degrees-of -freedom. 

This  report  presents  a  summary  of  the  theoretical  steps  involved 
in  the  development  of  the  element.  The  accuracy  of  the  element  is 
illustrated  by  the  solution  of  several  standard  problems  and  a 
comparison  of  results  with  other  thin  shell  elements.  The  FORTRAN 
77  listing  of  the  subroutines  which  form  the  basic  element 
matrices  contain  less  than  300  statements  and  is  presented  in 
order  to  illustrate  that  the  computer  implementation  of  the 
element  is  relatively  simple. 
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INTRODUCTION 

General  Background 

The  use  of  composite  materials  allows  for  the  efficient  design  of 
many  different  types  of  structural  systems.  One  of  the  major 
advantages  of  the  material  is  that  different  stiffness  and 
strength  properties  can  be  obtained  in  different  directions. 
Therefore,  more  efficient  structures  can  be  obtained  since  the 
material  can  be  concentrated  in  the  directions  of  maximum 
stresses. 

Most  of  the  existing  finite  element  programs  do  not  have 
sufficient  generality  to  consider  such  material  properties.  Also, 
in  the  case  of  thin  shell  structures  very  few  programs  have  the 
ability  to  consider  shells  in  which  the  bending  and  membrane 
forces  are  coupled.  In  addition,  problems  associated  with  the 
modelling  of  complex  shell  structures  with  thin  shell  elements 
exist  since  the  classical  thin  shell  formulation  does  not  have 
stiffness  terms  associated  with  the  normal  rotational  degrees  of 
freedom.  Therefore,  the  user  of  the  program  is  often  required  to 
add  artificial  members  to  a  finite  element  model  in  order  to  avoid 
numerical  instability  in  the  solution  of  the  finite  element 
system.  The  purpose  of  this  report  is  to  present  a  new  thin  shell 
element  which  is  sufficiently  robust  to  solve  the  above  mentioned 
problems . 

Recent  Research 

Within  the  past  two  years  several  papers  have  presented  methods 
which  introduce  a  normal  rotation  in  order  to  improve  the  membrane 
behavior  of  plane  elements.  Carpenter,  Stolarskl  and  Belytschko 
present  a  flat  triangular  shell  element  with  improved  membrane 
interpolation.  [1]  They  introduced  normal  mid-side  displacements 
of  the  constant  strain  triangle.  The  normal  displacements  are 
eliminated  and  node  rotations  are  introduced  by  the  use  of  a  cubic 
constraint  function  along  each  side  of  the  triangle.  A  one  point 
Integration  method  is  used  in  order  to  eliminate  membrane  locking 
within  the  elements.  The  element  yields  very  accurate 
displacements;  however,  the  element  is  rank  deficient  and  is 
unstable  for  certain  geometries. 
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Taylor  and  Slmo  applied  the  same  basic  approach  as  presented  in 
reference  [1]  to  improve  the  membrane  behavior  of  quadrilateral 
elements  [2] .  For  many  problems  excellent  displacements  and 
stresses  are  obtained.  However,  for  shell  structures  such  as  a 
twisted  beam  the  displacements  become  very  large  as  the  mesh  was 
refined.  In  addition,  the  flat  quadrilateral  element  does  not 
accurately  model  many  common  types  of  shell  geometries.  Also,  the 
DKQ  formulation  was  used  to  form  the  bending  stiffness  which  has 
proven  to  be  not  as  accurate  as  the  DKT  formulation. 

Bergan  and  Felippa  have  developed  a  triangular  membrane  element 
with  normal  rotational  degrees-of-freedom  [3] .  The  formulation  is 
based  on  the  "free  formulation"  and  uses  the  continuum-mechanics 
definition  of  rotation.  The  element  passes  the  patch  test  and  is 
stable  for  all  applications.  The  element  produces  good  values  of 
displacements;  however,  the  values  for  stresses  are  poor  compared 
to  the  values  obtained  from  the  Taylor  quadrilateral  [4] . 

The  three  elements  previously  mentioned  have  not  been  used  for 
thin  shells  in  which  the  membrane  and  bending  forces  are  coupled. 
Therefore,  one  of  the  purposes  of  this  report  is  to  develop  an 
element  which  has  a  consistent  formulation  for  both  the  bending 
and  membrane  behavior.  In  addition,  the  problems  with  the 
instability  associated  with  normal  rotational  degree-of-f reedom 
will  be  studied  and  a  simple  technique  is  suggested  in  order  to 
avoid  this  problem. 
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BASIC  BQUATI(»fS  -  ORTBOTROPIC  MATERIALS 

The  18  X  18  triangular  shell  element  stiffness  matrix  for  a  stiffened 
composite  material  as  shown  in  figure  1  can  be  directly  calculated 
from  the  following  well-known  equation: 

E  ■  I  IT  B  B  dA  (1) 

The  6x6  O  matrix  relates  the  forces  to  the  deformations  which  are 
associated  with  a  differential  element  of  area  dA.  Including  thermal 
deformations  the  force-deformation  relationship  can  be  expressed  by 
the  following  matrix  equation: 

1  -  D  e  +  f .  (2) 

where 

f  ^  [  mi  1  ms  t  mi  t  f  i  t  f  1 1  f  1  s  ]  ( 3 ) 

and 

A*  *  [  kii  kti  kit  cti  ttt  tis  ]  (4) 

The  positive  definition  of  these  forces  and  deformations  is 
illustrated  in  figure  2. 


Normally  the  matrix  0  cannot  be  defined  directly  for  complex 
materials.  However,  the  inverse  can  normally  be  easily  calculated 
from  the  basic  principles  of  mechanics  or  determined  experimentally. 
Therefore,  the  terms  for  are  normally  specified  as  input  to  a 

computer  program.  The  numerical  values  of  D  are  then  evaluated  within 
the  element  stiffness  subroutine.  Hence,  the  basic  behavior  of  the 
thin  shell,  including  thermal  deformations,  is  expressed  in  the 
following  form: 
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Or  in  terms  of  matrix  notation: 

e.  -  D-»  f  +  <•  (5) 

Where  dT  is  the  temperature  change  and  at  to  cu  are  the  measured 
thermal  expansion  coefficients.  Hence,  the  thermal  forces  indicated 
in  equation  (2)  are  calculated  from: 

f«  "  -  D  (6) 

Each  flexibility  term  in  equation  (6)  has  a  direct  physical  meaning. 
For  example,  the  term  Cia  is  the  curvature  kn  due  to  a  unit  force, 
fat  »  1.  Whereas,  the  term  Cat  is  the  strain  C]^]^  at  the  reference 
plane  caused  by  the  application  of  a  unit  bending  moment,  mi i  ■  1. 
It  is  apparent  that  these  terms  can  be  best  determined  experimentally 
for  complex  composite  materials.  Also,  the  values  of  these  terms  are 
dependent  on  the  definition  of  the  reference  plane  which  must  be 
defined  at  the  same  time  as  the  flexibility  terms  are  determined. 
For  the  special  case  of  constant  thickness  isotropic  shells  the 
mid-surface  is  the  logical  reference  plane  and  the  terms  Cij  and 
several  other  flexibility  terms  are  zero  in  equation  (5) . 

In  equation  (1)  the  6x18  B  matrix  defines  the  relationship  between 
the  deformation  terms  and  the  node  displacements  ▼  in  the  local 
1,2,3  system.  Or,  in  matrix  form: 

1.  -  B  ▼  (7) 

which  can  be  written  in  submatrix  form  as 

▼  (8a) 

•  Bp  ▼  (8b) 

where  the  "p”  and  "m"  indicate  the  plate-bending  and  membrane  terms 

respectively. 
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BENDING  APPROXIMATION  -  THE  DRT  ELEMENT 

The  development  of  the  B»  matrix  la  baaed  on  the  atandard  DAT  element 
[5] .  Becauae  the  bending  and  membrane  behavior  are  coupled  the  DXT 
formulation  will  be  aummarized  here  in  order  that  the  B*  matrix  will 
be  developed  with  conaiatent  approximationa  in  the  next  section  of 
this  report. 

The  DKT  element  is  based  on  the  independent  expansion  of  the  inplane 
rotations  of  the  reference  surface  for  a  6  node  triangle  which  is 
shown  in  figure  3.  If  the  local  1  and  2  directions  are  indicated  by 
the  local  x  and  y  coordinates  the  finite  element  approximation  is 
written  in  the  following  form: 

Ps  ■  ai  +  atx  +  aty  ♦  Ka4X*  +  asxy  Xaty* 

(9) 

Pr  "a?  +  a«x  +  a»y  +  Kat*x*  +  aiixy  +  Xaity* 

The  six  constants  at  to  a« ,  as ,  can  be  expressed  in  terms  of  the  six 
node  rotations  Pxt  to  Pxt  by  an  inversion  of  a  6x6  matrix  which  will 
produce  an  equation  of  the  following  form: 


1*  ■ 

1  fts 

(10) 

The 

same 

6x6  matrix,  H,  will  relate  the  constants  st  to  ait. 

o 

m 

the 

node 

rotations  fir . 

From  the  theory  of  thin  plates  the  curvature-displacement 
relationships  are  defined  by  the  following  equations: 


kss  •  Ps  ,1  ■  as  a4X  asy  (11a) 

ksr  ■  Pr»F  *  as  >  aitx  atsy  (lib) 

k*f  ■  Ps  fF  +  Pf  f* 

.s 

*  as  asx  *  asy  as  -»>  aisx  *  any  (11c) 

Or,  in  the  following  matrix  form: 

As  *  Gf  a  ( 12 ) 


(a)  BASIC  ROTATIONAL  UNKNOWNS 


Ljj  is  the  length  from  I  to  J 


I  M  J 


(b)  DISPLACEMENTS  ALONG  TYPICAL  SIDE 


Figure  3.  DISPLACEMENT  APPROXIMATION  -  DRT  ELEMENT 


Where 


“fa*  +a«  +  a»  •♦■a4  ■♦•aa  •♦■aa  + 


and 


& 


aa  +  aa  +  ai  a  +  ai  i  +  at  a  ] 


lOxyOOOOOO 

OOOOOOlOxy 

OlOxylOxyO 


(13) 


(14) 


The  six  rotational  degrees-of-freedom  which  are  associated  with  a 
typical  side  I-J  of  the  element  are  shown  in  figure  (3a) .  These 
rotations  can  be  transformed  to  a  n-s  reference  system  which  is 
parallel  and  normal  to  a  typical  element  side  as  shown  in  figure  (3b). 
The  basic  DRT  constraints  are  enforced  as  follows: 


1.  The  mid-side  rotation  B»m  is  set  to  the  average  of  the  values 
at  point  I  and  J.  Or, 

e»m  •  {  9ml  *  Sms  )  /  2  (15) 

2.  The  s-displacements,  above  and  below  the  reference  plane,  and 
the  normal  displacement  in  the  z-direction  are  forced  to  be 
cubic  functions  since  the  transverse  shear  strain  is  in  the 
s-direction  is  set  to  zero.  Therefore,  the  mid-side  normal 
rotation  must  satisfy  the  following  equation: 

eiUB  ■  -  {9mi  *  9ml)/4  -  3{wi  Wj)/(2Lij)  (16) 

The  constraint  specified  by  equation  (15)  will  force  the  normal 
displacement  along  the  element  sides,  above  and  below  the  reference 
plane,  to  be  linear  function.  Hence,  displacement  and  slope 
compatibility  is  satisfied  along  the  sides  of  all  elements.  Since  no 
attempt  is  made  to  set  the  transverse  shear  strains  to  zero  within  the 
element  the  name  Discrete  Klchoff  Triangle  was  selected  as  the  name  of 
this  element. 


9 


I 


Equation  (IS)  and  (16)  can  be  aununarized  in  matrix  form  as 


0.  ' 

'l/2  o' 

Os  I 
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■  1/2 

o' 

Os  j 

+  1 . 5/Li  J 
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'Wi  ■ 

Oh 

0  -1/4 

On  I 
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-1/4 
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-1 

1 

Wj 

The  relationship  between  the  N-S  and  X-Y  coordinate  systems  are 


Os 

Cos5 

SinS 
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-SinS 
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Or 
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Ox 
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On 

Equation  (17)  can  now  be  written  in  the  X-Y  system  as 


Ox 

Tl 

T2  ' 

'  Oxi  ' 

+ 
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Wj 

Where 


(17) 


(18a) 


(18b) 


(19) 


T1  »  0.5  Cos* 5  -  0.25  Sin* 5 
T2  ■  1.5  SinS  Cos5 
T3  -  0.5  Sin* 5  -  0.25  Cos* 6 
Ts  •  1.5  Li j  SinS 
Tc  «  1.5  Li /  Cos5 

The  six  rotational  deqrees-o£-freedom  associated  with  the  three 
mid-side  nodes  of  the  triangle  can  be  eliminated  by  the  application  of 
equation  (19) .  These  transformations  can  be  summarized  by  a  matrix 
equation  of  the  following  form: 


e  -  Tr  w 


(20) 


where  Tr  is  a  12x9  matrix  and  W  represents  a  9x1  vector  of  Ox i ,  9v  i 
and  wi  for  the  three  nodes  of  the  triangle. 
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MEMBRANE  APPROXIMATION 

The  membrane  behavior  of  the  triangular  shell  element  is  based  on  the 
basic  six  node  quadratic  membrane.  If  the  12  coefficents  are  defined 
by  the  12  x  1  vector  b  the  membrane  strains  can  be  expressed  in  the 
following  form: 

JL"  “  2"  fe  (21) 

As  in  the  case  of  the  DKT  plate  bending  element  the  three  mid-side 
displacements  are  rotated  to  the  local  N-S  coordinate  of  each  side  as 
shown  in  figure  (4) .  In  order  to  maintain  displacement  compatibility 
between  element  the  displacement  us  is  assumed  to  be  linear  along  each 
side  and  the  displacement  uh  is  a  cubic  function.  These  assumptions 
can  be  summarized  by  the  following  equations  for  the  displacements  at 
the  mid-side  nodes: 

Us  ■  (xis  1  +  lu  j  )  /2 

(22) 

Un  =  (umi  +  Ui»j)/2  +  Li  j  (Bz  I  -  Bzj)/8 

These  equations  can  be  written  in  terms  of  the  X-Y  coordinate  system 
as 
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The  six  translational  displacements  at  the  midside  nodes  can  now  be 
eliminated  and  three  rotational  unknowns  are  introduced  at  the 
vertices  by  the  direct  application  of  equation  (23) .  The  same  basic 
approach  which  was  used  in  the  DKT  formulation  is  now  applied  in  order 
to  form  the  matrix  equation  of  the  following  form: 

b  s  Tn  tt  (24) 

Therefore,  the  three  membrane  strains  can  be  written  in  terms  of  the 
nine  node  displacements  as 

Cm  ■  Gn  ^  tt  (25) 


It  is  now  possible  to  evaluate  the  24  x  24  element  stiffness  by  the 
direct  application  of  equation  (1) . 


COMPUTBR  PROGRAM  IMPLEMENTATION 


The  18  X  18  triangular  element  stiffness  matrix,  for  the  element  shown 
in  figure  5,  is  given  by  equation  (1).  Where  the  the 

strain-displacement  matrix  B  can  now  be  written  in  terms  of  the 
bending  and  membrane  submatrices,  or 


Gp  0  Tp 

0 

B  - 

■  G(x,y)  T 

(26) 

0  Gm  0 

Tm 

Since  the  T  matrix 

is 

not  a  function  of  x  and  y  it 

is  possible  to 

rewrite  equation  (1) 

in 

the  following  form: 

K  «  TT  [  G»  D  G  dA 

T 

(27) 

The  matrix  G  is  very  sparse  and  contains  only  the  terms  1,  x  and  y; 
therefore  the  integral  cab  be  evaluated  directly  with  a  minimum  of 
numerical  effort  as  Illustrated  by  the  FORTRAN  listing  given  in 
Appedix  A.  In  addition,  an  integral  reduction  factor  can  be  used  on 
selective  terms  in  order  to  improve  the  membrane  performance  as 
suggested  in  reference  (3) . 

NUMERICAL  EXAMPLES 


In  order  to  illustrate  the  behavior  of  the  element  and  to  compare  the 
results  with  other  shell  elements  several  examples  will  be  presented. 

Cantilever  Beam 

The  beam  shown  in  figure  6  is  idealized  by  a  1x4  rectangular  mesh  and 
is  subjected  to  a  load  of  40  kips  at  the  tip  of  the  cantilever.  The 
theoretical  displacement  at  the  tip,  including  shearing  deformation, 
is  0.3558  inches. 
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The  Taylor  quadrilateral  element  shell  yields  a  displacement  of  0.3467 
inches;  or  an  error  of  -1.02  percent.  Note  that  the  rotations  at  the 
base  of  the  cantilever  are  set  to  zero  which  is  inconsistent  with  the 
existence  of  shearing  deformations. 

The  element  presented  In  this  report,  TSHELL,  was  used  to  model  this 
beam  with  two  triangles  to  form  each  quadrilateral.  The  completely 
integrated  element  produces  a  displacement  of  0.2695  inches,  or  an 
error  of  -24.3  percent.  With  a  reduced  integration  factor  of  0.5  the 
tip  displacement  is  0.3726,  or  an  error  of  '»-4.5  percent. 

For  all  example  problems  presented  in  this  report  a  reduced 
integration  factor  of  0.5  is  used.  The  use  of  the  reduced  integration 
factor  has  the  major  advantage  over  one  point  integration  is  that  a 
"rank  deficiency"  is  not  introduced  into  the  element. 


Figure  6.  CANTILEVER  BEAM  EXAMPLE 


14 


I 

Twisted  Bea» 

The  twisted  beam  shown  in  figure  7  has  become  a  standard  test  problem 
I  for  thin  shell  elements  [6] .  Table  1  summarizes  the  results  obtained 

using  four  different  elements.  It  is  important  to  note  that  the 
SHELL  element  does  not  converge  as  the  mesh  is  refined.  The  reason 
for  this  unacceptable  behavior  is  because  the  element  is  flat  and  it 
I  cannot  model  the  twisted  surface  accurately.  The  new  triangular 

element  gives  good  results  for  this  problem. 


WIDTH*  1. 1 
DEPTH  *0.32 
TWIST  «  90* 
E«29.0x  10® 
POISSON  RATIO  *0.22 
MESH  •  12  X  2 
Unit  loads  of  end 

Figure  7.  TWISTED  BEAM  EXAMPLE 


Table  1.  Deflection  Under  Load  For  Twisted  Beam 


Element  Type  and  Mesh 

IN-PLANE 

error 

OUT-OF- PLANE 

error 

Reference  Values 

0.005424 

0 

0.001754 

0 

NASTRAN  QUAD4 

0.005386 

-0.7 

0.001727 

-1.5 

NASTRAN  QUADS 

0.005413 

-0.2 

0.001750 

-0.2 

SHELL  (2x12) 

0.005779 

*6.  A 

0.001993 

+13.7 

"  (4x24) 

0.006849 

426.2 

0.002750 

+56.9 

TSHELL  (2x12) 

0.005390 

-0.6 

0.001717 

-2.1 

**  (4x24) 

0.005399 

-0.5 

0.001735 

-1.1 
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Scord«ll0-Lo  Roof 

Th*  shell  structure  shown  In  figure  8  is  a  standard  test  problem  [6] . 
Table  2  summarizes  the  maximum  displacement  obtained  using  different 
elements  and  meshes.  For  this  problem  the  SHELL  element  results  are 
very  good.  However,  the  TSHELL  results  appear  to  converge  very 
slowly.  However,  for  even  the  coarse  mesh,  the  results  are  of 
acceptable  accuracy  for  normal  engineering  analysis. 


THICKNESS  «0.25 
E«  4.32x10® 

POISSON  RATIO*  0.0 
loading :  90/unif  ores  in  z  dir. 
u  *  w  *  0  on  curved  edge 
MESH  NxN  onquodrant 


Figure  8.  SCOROELIS-LO  CYLINDRICAL  SHELL 


Table  2.  Maximum  Displacement  Of  Scordelis-Lo  Roof 


Element 

Type  and  Mesh 

VALUE 

ERROR 

Reference  Value 

0.3086  ft. 

0.0  % 

SHELL 

N-4 

•»-5.2  ! 

N>6 

-♦•1.7 

N-8 

■••0.6  ! 

N-10 

•fO.2  ! 

TSHELL 

N-4 

0.3036 

-1.6 

N-8 

0.2982 

-3.4 

N-12 


0.2997 


-2.9  % 
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FINAL  REMARKS 

A  new  anisotropic  triangular  shell  element,  with  normal  rotations  at 
the  nodes  has  been  developed.  The  element  can  be  connected  directly 
to  nodes  with  beam  elements  without  special  consideration. 

The  general  accuracy  of  the  element  has  been  demonstrated.  Care 
should  be  taken  if  the  element  is  used  to  model  shells  which  have 
double  curvature. 
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APPENDIX  A  -  FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 


C - SHELLT 

SUBROUTINE  SHELLT (S,T,F.D,C,V,X,Y, AREA , PRES , TEMP ) 
IMPLICIT  REAL*8  (A-H,0-Z) 

C -  INFORMATION  CALCULATED  - 

C  S  -  18  X  18  STIFFNESS  MATRIX  IN  X-Y-Z  SYSTEM 

C  T  -18x18  FORCE-TRANSFORMATION  MATRIX 

C  F  «  18  X  2  FORCES  DUE  TO  PRESSURE  AND  TEMP. 

C -  INFORMATION  SPECIFIED  - 

C  D  «  6  X  6  MATERIAL  PROPERTY  MATRIX 

C  LOCAL  FORCES  ARE  ASSUMMED  TO  BE  IN  THE 

C  ORDER  Mil,  M22,  M12,  Nil,  N22  and  N12 

C  C  -  6  X  1  MATRIX  OF  THERMAL  EXPANSION  TERMS 

C  V  -  THE  DIRECTION  COSINE  ARRAY  WHICH  RELATES  THE 

C  LOCAL  1-2-3  TO  THE  GLOBAL  X-Y-Z  SYSTEM 

C  X,  Y,  Z  GLOBAL  NODE  COORDINATES 

C  AREA  -  ELEMENT  AREA 

C  PRES  «  AVERAGE  SURFACE  PRESSURE 

C  TEMP  -  AVERAGE  TEMPERATURE  CHANGE 

C - 

C  WRITTEN  BY  EDWARD  L.  WILSON,  JAN.  -  JUNE  1987 

C - 

DIMENSION  V(4,4) ,RL(24,3) ,XY(3) ,E(6) ,H(6,6) , 


.  S  (18, 18) ,T (18,18) ,D (6,6) ,X(6) ,Y(6) ,AP (10,12) , 

.  AM(10,12) ,A(20,18) , AIN (3, 3) ,G (20 , 20) , JP (9) , 

.  JM(9) ,TT(20,18) ,F(18,2) ,C(6) 

DATA 

.  JP  /4, 5, 10, 11, 16, 17, 3, 9, 15/, 

.  JM  /I, 2, 7, 8, 13, 14, 6, 12, 18/ 

DATA  KL  / 

.  1,1,1,  2,2,2,  3, 3, 3, 3, 3, 3,  4,4,4,  5,5,5,  6, 6, 6, 6, 6, 6, 

.  1,3,4,  7,9,10,  2, 4, 5, 6, 8, 9,  11,13,14,  17,19,20,12,14,15,16,18,19, 
.  1,2,3,  1,2,3,  1,2, 3, 1,2, 3,  1,2,3,  1,2,3,  1,2, 3, 1,2, 3/ 

NT  »  24 

C -  CHANGE  TO  LOCAL  COORDINATES  SYSTEM  - 

XC  -  (  X(l)  +  X(2)  +  X(3)  )  /  3. 

YC  »  (  Y(l)  +  Y(2)  +  Y(3)  )  /  3. 

DO  20  1-1,3 
X(I)  -  X(I)  -  XC 

Y(I)  -  Y(I)  -  YC 

DO  20  J-1,3 
20  AIN(I,J)  -  0.0 

C -  COMPUTE  INTEGRALS  - 

RED  -  0.50* AREA 
AIN(1,1)  -  AREA 

AIN(2,2)  — RED*(  X(1)*X(2)  ^  X(2)*X(3)  X(3)*X(1)  )  /  6. 

AIN(2,3)  -  RED*(  X(1)*Y(1)  X(2)*Y(2)  +  X(3)*Y(3)  )  /  12. 

AIN(3,3)  — RED*(  Y(1)*Y(2)  Y(2)*Y(3)  +  Y(3)*Y(1)  )  /  6. 

AIN(3,2)  -  AIN(2,3) 


-  19  - 

yORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C -  FORM  “H"  ARRAY  FOR  GENERAL  6-NODE  TRIANGLE  - 

CALL  FORMH(H,A,X,Y) 

C -  FORM  "AP"  ARRAY  FOR  PLATE  ELEMENT  - 

CALL  FORMAP(AP,H,X,Y) 

C -  FORM  "AM"  ARRAY  FOR  PLANE  ELEMENT  - 

CALL  FORMAM(AM,H,X,Y) 

C -  FORM  20x20  "G"  ARRAY  - 

DO  160  K«l,20 
DO  160  L«1,K 
160  G(K,L)  «  0.0 
C 

DO  200  N«1,NT 

I  «  KL(N,1} 

K  »  KL(N,2) 

II  -  KL(N,3) 

DO  190  M«1,NT 
J  -  RL(M,1) 

IF  (D(I. J) .EQ.0.0)  GO  TO  190 
L  «  RL(M,2) 

IF  (L.GT.K)  GO  TO  190 
JJ  >  KL(M,3) 

G(K,L)  -  G(K.L)  +  D(I,J)*AIN(II,JJ) 

190  CONTINUE 
200  CONTINUE 
C 

DO  210  K-1.20 
DO  210  L»1,K 
210  G(L.K)  «  G(K.L) 

C -  FORM  20x18  "A"  ARRAY  - 

DO  300  J-1,18 
F(J,1)  -  0.0 
F(J.2)  -  0.0 
DO  300  1-1,20 
300  A(I,J)  -  0.0 
C 

DO  305  J-1,9 
JJ  -  JP(J) 

DO  305  1-1,10 
305  A(I,JJ)  -  AP(I,J) 

C 

DO  310  J-1,9 
JJ  -  JM(J) 

DO  310  1-1,10 
310  A{H'10,JJ)  -  AM(I,J) 

C -  ROTATE  TO  GLOBAL  X,  Y,  Z  SYSTEM  - 

DO  350  N-1,6 
NZ  -  3*N 
NY  -  NZ  -  1 
NX  -  NY  -  1 
DO  350  1-1,20 

XX  -  A(I,NX)»V(1,1)  +  A(I,NY) *V(1,2)  +  A(I,NZ) *V(1,3) 

YY  -  A(I,NX) •V(2,l)  +  A(I,NY)*V(2,2)  +  A(I,NZ) *V(2,3) 

ZZ  -  A(I,NX) *V(3,1)  +  A(I,NY)»V(3,2)  +  A(I,NZ) »V(3, 3) 

A(I,NX)  -  XX 
Ad, NY)  -  YY 
350  Ad,NZ)  -  ZZ 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C -  FORM  18x18  ELEMENT  STIFFNESS  - 

DO  380  I«l,20 
DO  380  J=l,18 

CALL  DOTP(G(l,I)  ,A(1,  J)  ,SX»1,20) 

380  TT(I,J)  »  SUM 
C 

DO  400  I«l,18 
DO  400  J«1,I 

CALL  DOTP{A(l,I) ,TT{1, J) ,SUM,20) 

S(1,J)  »  SUM 
400  S(J,I)  -  SUM 

C -  FORM  FORCE-DISPLACEMENT  TRANSFORMATION  ARRAY 

NO  «  0 

DO  500  N»l,3 
XY(1)  «  1.0 
XY(2)  =  X(N) 

XY(3)  «  Y(N) 

C 

DO  450  K«l,18 
C 

DO  410  I«l,6 
410  E(I)  «  0.0 
C 

DO  420  M»1,NT 

I  »  KL(M,1) 

L  «  KL(M,2) 

J  »  KL(M,3) 

420  E(I)  a  E(I)  +  XY{J)*A(L,K) 

C 

DO  440  1-1,6 
SUM  -  0.0 
DO  430  J«l,6 

430  SUM  -  SUM  +  D(J,I)*E(J) 

II  »  NO  +  I 
440  T(II,K)  -  SUM 

C 

450  CONTINUE 
C 

500  NO  -  NO  6 

C -  FORM  THERMAL  FORCES  - 

IF  (TEMP. EQ. 0.0)  GO  TO  625 
DO  600  1-1,6 

CALL  DOTP(D(l,I) ,C(1) ,E(I) ,6) 

600  CONTINUE 

DO  610  1-1,6 
610  C(I)  »  -  TEMP*E(I) 

C 

DO  620  1-1,18 
F(I,2)  -  C(1)*A(1,I) 

F(I,2)  -  C(2)*A(7,I) 

F(I,2)  -  C(3)*(  A(2,I)  A(6,I)  ) 

F(l,2)  -  C(4)*A(11,I) 

F(I,2)  -  C{5)*A(17,I) 

620  F(I,2)  -  C(6)»(  A(12,I)  +  A(16,I)  ) 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C -  CALCULATE  PRESSURE  FORCES  - 

625  IF ( PRES. EQ. 0.0}  GO  TO  800 
FORCE  «  -  AREA*PRES/3.0 
DO  630  1-1,18,6 
F(I  ,1)  -  PORCE*V(l,3) 

F(  1+1,1)  »  FORCE* V( 2, 3) 

630  F(I+2,1)  »  FORCE*V(3,3) 

C 

800  RETURN 

- - 

END 

C  FORMH 

SUBROUTINE  FORMH (H , B , X,Y) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  H(6,6),B(6.6).Y{4) ,X(4) 

C -  FORM  COEFFICIENT  MATRIX  FOR  6  NODE  TRIANGLE  - 

X(4)  -  (  X(l)  +  X(2)  )  /  2. 

X(5)  »  (  X(2)  +  X(3)  )  /  2. 

X(6)  -  (  X(3)  +  X(l)  )  /  2. 

Y(4)  -  (  Y(l)  +  Y(2)  )  /  2. 

Y(5)  »  (  Y(2)  +  Y(3)  )  /  2. 

Y(6)  -  (  Y(3)  +  Y(l)  )  /  2. 

C 

DO  100  1-1,6 
H(1,I)  -  1.0 

H(2,I)  -  X(I) 

H(3,I)  -  Y(I) 

H(4,I)  «  X(I)*X{I)  /  2. 

H(5,I)  -  X(I)*Y(I) 

100  H(6,I)  «  Y(I)*Y(I)  /  2. 

C -  INVERT  TO  FORM  H  MATRIX  - 

DO  200  1-1,6 

DO  200  J-1,I 

SUM  -  0.0 
DO  190  K-1,6 

190  SUM  -  SUM  +  H(I,K) *H{J,K) 

B(J,I)  -  SUM 
200  B(I,J)  -  SUM 
C 

CALL  SYMSOL(B,H,6,6,0) 

C 

RETURN 

END 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C - FORMAP 

SUBROUTINE  FORMAP (AP , H , X , Y) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  AP(10,12) ,H(6.6) ,Y(4) .X(4) ,IT(4) 

DATA  IT  /I. 3, 5,1/ 

C -  FORM  12  DOF  MATRIX  - 

DO  240  I«2,6 
K  =  0 

DO  240  Jsl,6 
K  »  K  +  1 
AP(I-1,K)  »  0.0 
AP(I+4,K)  =  H{I,J) 

K  »  K  +  1 

AP(I-1,K)  =  -  H{I,J) 

240  AP(I+4,K)  »  0.0 

C -  ELIMINATION  OF  4,5,6  MID-SIDE  ROTATIONS  - 

X(4)  =  X(l) 

Y(4)  -  Y(l) 

DO  300  N-1,3 

DX  =  X(N+1)  -  X(N) 

DY  =  Y(N+1)  -  Y(N) 

XL  -  DSQRT(  DX*DX  +  DY*DY  ) 

S  »  DY  /  XL 
C  «  DX  /  XL 
T1  «  C*C/2.  -  S*S/4. 

T2  =  0.75*S*C 

T3  =  S*S/2.  -  C*C/4. 

XX  =  1.5/XL 
TC  «  XX*C 
TS  «  XX*S 
C 

11  »  IT(N) 

12  «  II  +  1 
J1  »  IT(N+1) 

J2  *  J1  +  1 

L2  »  6  +  N  +  N 
LI  -  L2  -  1 
I  -  N  ^  6 
J  -  I  +  1 
IF  (N.EQ.3)  J  -  7 
C 

DO  300  K»l,10 

AP(K,I1)  »  AP(K,I1)  +  AP{K,L1)*T1  +  AP(K,L2)*T2 

AP(K,I2)  »  AP(K,I2)  4-  AP{K,L1)*T2  +  AP{K,L2)*T3 

AP(K,J1)  «  AP(K,J1)  +  AP{K,L1)*T1  +  AP(K,L2)*T2 

AP(K,J2)  »  AP(K,J2)  +  AP(K,L1)*T2  +  AP{K,L2)»T3 

W  »  -  AP(K,L1)*TS  ♦  AP(K,L2)*TC 
AP(K,L1)  -  0.0 

AP(K,L2)  -  0.0 

AP(K,I)  »  AP(K,I)  +  W 
AP(K,J)  «  AP(K,J)  -  W 
300  CONTINUE 
C 

RETURN 

END 


FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C - FORMAM 

SUBROUTINE  FORMAM (AM. H.X.Y) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

DIMENSION  AM(10,12) ,H(6,6) ,Y(4) ,X(4) ,IT(4) 

DATA  IT  /I, 3, 5,1/ 

C -  FORM  12  DOF  MATRIX  - 

DO  240  Is2,6 
K  »  0 

DO  240  Jsl.e 
K  »  K  +  1 
AM(I-1,K)  =  H(l.J) 

AM(I-«-4,K)  -  0.0 
X  «  K  +  1 
AM(I-1,K)  »  0.0 

240  AM(I+4,K)  »  Hd.J) 

C -  ROTATE  NODE  4,  5,  6  MID-SIDE  DISPLACEMENTS  - 

X(4)  =  X(l) 

Y(4)  -  Y(l) 

C 

DO  300  N=l,3 

DX  =  X{N+1)  -  X(N) 

DY  -  Y(N+1)  -  Y(N) 

C 

NY  «  2*N  +  6 
NX  =  NY  -  1 
IX  *  IT(N) 
lY  =  IX  +  1 
JX  ■  IT(N+1) 

JY  «  JX  +  1 
I  «  N  +  6 
J  -  I  +  1 
IF  (N.EQ.3)  J  ■  7 

C -  ELIMINATE  MID-SIDE  DISPLACEMENTS  - 

DO  260  K»l,10 

TT  =  .125*(  -  AM(K,NX)*DY  AM(K,NY)*DX  ) 

AM(K,IX)  «  AM(K,IX)  AM(K,NX)/2. 

AMIK.IY)  »  AM(K,IY)  +  AM{K,NY)/2. 

AM(K,JX)  «  AM(R,JX)  *  AM(K,NX)/2. 

AM(K,JY)  «  AM(K,JY)  +  AM(K,NY)/2. 

AM(K,NX)  -  0.0 
AM(K.NY)  -  0.0 
AM(K,I)  «  AM(K,I)  +  TT 
AM(K,J)  -  AM(K,J)  -  TT 
260  CONTINUE 
C 

300  CONTINUE 
C 

RETURN 

END 


o  o  o 
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FORTRAN  LISTING  OF  TRIANGULAR  SHELL  ELEMENT 

C - LOCALT 

SUBROUTINE  LOCALT ( XYZ , X , Y , V , AREA , I AXI S ) 

IMPLICIT  REAL* 8  (A-H,0-Z) 

DIMENSION  XYZ(3,3) ,X(6) ,Y(6) ,V(4,4) 

-  TRANSFORM  TO  LOCAL  COORDINATE  SYSTEM  - 

WRITE  (*,3000)  XYZ 
3000  FORMAT  (3F15.4} 

IF  (lAXIS.NE.O)  THEN 
DO  100  I»l,3 
100  V(I,4)  »  0.0 

I  -  lABS(IAXIS) 

V(I,4)  *  lAXIS/I 
END  IF 

C -  DEFINE  LOCAL  1,2,3  REFERANCE  SYSTEM - 

CALL  VECTOR (V (1,1) , XYZ (1,1) ,XYZ(1,2) , XYZ (1,3) , 

XYZ(2,1) ,XYZ(2,2) .XYZ(2,3) ) 

CALL  VECTOR(V(l,2) ,XYZ(1,1) ,XYZ(1,2) ,XYZ(1,3) , 

XYZ{3,1) ,XYZ(3,2) ,XYZ(3,3) ) 

CALL  CROSS(V(l,l) ,V(1,2) ,V(1,3) ) 

IF  (lAXIS.NE.O)  CALL  CROSS(V(l,4) ,V(1,3) ,V(1,1) ) 

CALL  CROSS(V(l,3) ,V(1,1) ,V(1,2) ) 

C 

DO  5  N=l,3 

X(N)  »  XYZ(N,1)  *V(1,1)  XYZ(N,2)  *V(2,1)  +  XYZ  (N,  3)  *V(3 , 1) 
5  Y(N)  «  XYZ(N,1)*V(1,2)  +  XYZ (N, 2 ) *V ( 2 , 2 )  +  XYZ (N, 3 ) *V ( 3 , 2) 

C -  CALCULATE  AREA  OF  ELEMENT  - 

AREA  «  (  X(2)*Y(3)  -  X(3)*Y(2)  +  X(3)*Y(1) 

-  X(1)*Y(3)  +  X(1)*Y(2)  -  X(2)*Y(1)  )  /  2.0 
IF  ( AREA. LE. 0.0)  THEN 
PAUSE  '  ZERO  OR  NEGATIVE  AREA  ' 

RETURN 
END  IF 

C - 


RETURN 

END 
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